function x = erfcinvO( y )
%erfcinvO: Computes the inverse of complementary error function. 

% erfcinvO: A different function from Matlab's erfcinv. Here 'O' is used to
% differentiate these two functions

if y < 0 || y >2
    error ('Inverse complementary function not defined outside [0, 2]');
elseif y == 0
    x = inf;
elseif y == 2
    x = -inf;
else
    if y >= 0.0485 && y <= 1.9515
        q = y - 1;
        r = q * q;
        x = (((((0.01370600482778535*r - 0.3051415712357203)*r + 1.524304069216834)*r - 3.057303267970988)*r + 2.710410832036097)*r - 0.8862269264526915) * q /...
          (((((-0.05319931523264068*r + 0.6311946752267222)*r - 2.432796560310728)*r + 4.175081992982483)*r - 3.320170388221430)*r + 1.0);
    elseif y < 0.0485
        q = sqrt (-2.0 * log (y / 2.0));
        x = (((((0.005504751339936943*q + 0.2279687217114118)*q + 1.697592457770869)*q + 1.802933168781950)*q + -3.093354679843504)*q - 2.077595676404383) / ...
          ((((0.007784695709041462*q + 0.3224671290700398)*q + 2.445134137142996)*q + 3.754408661907416)*q + 1.0);
    elseif y > 1.9515
        q = sqrt (-2.0 * log (1.0 - y / 2.0));
        x = (-(((((0.005504751339936943*q + 0.2279687217114118)*q + 1.697592457770869)*q + 1.802933168781950)*q + -3.093354679843504)*q - 2.077595676404383) / ...
                     ((((0.007784695709041462*q + 0.3224671290700398)*q + 2.445134137142996)*q + 3.754408661907416)*q + 1.0));
    else 
        x = 0;
    end
    u = (erfcO (x) - y) / (-2.0 / sqrt (pi) * exp (-x * x));
    x = x - u / (1 + x * u);   
end

end

